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Abstract 

We present a constitutive model capturing some of the experimen- 
tally observed features of soft biological tissues: nonlinear viscoelas- 
ticity, nonlinear elastic anisotropy, and nonlinear viscous anisotropy. 
For this model we derive the equation governing rectilinear shear mo- 
tion in the plane of the fiber reinforcement; it is a nonlinear partial 
differential equation for the shear strain. Specializing the equation 
to the quasistatic processes of creep and recovery, we find that usual 
(exponential-like) time growth and decay exist in general, but that for 
certain ranges of values for the material parameters and for the angle 
between the shearing direction and the fiber direction, some anoma- 
lous behaviors emerge. These include persistence of a non-zero strain 
in the recovery experiment, strain growth in recovery, strain decay in 
creep, disappearance of the solution after a finite time, and similar 
odd comportments. For the full dynamical equation of motion, we 
find kink (traveling wave) solutions which cannot reach their assigned 
asymptotic limit. 

1 Introduction 

Many biological, composite, and synthetic materials must be modeled as 
fiber-reinforced nonlinearly elastic solids. Hence, the anisotropy due to the 
presence of collagen fibers in many biological materials has been studied ex- 
tensively within the constitutive context of fiber-reinforced materials by sev- 
eral authors (see for example Humphrey (2002) and the references therein.) 
In nonlinear elasticity, the macroscopic response of an anisotropic material 
is given in terms of a strain-energy function, which itself depends on a set of 
independent deformation invariants. This formulation captures a great vari- 
ety of phenomena related to the behavior of fiber-reinforced materials such 
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as inter alia, the examination of fiber instabilities, using loss of ellipticity 
(see Merodio and Ogden 2002, 2003, and the references therein). 

Generally speaking, a reinforcement is added to a given material with 
the aim of avoiding a possible failure under operating conditions. Therefore 
it is important to develop a detailed study showing how to introduce rein- 
forcements in a material in order to control the possible development of a 
boundary layer structure. Our goal here is to provide a first step in this di- 
rection. We make several simplifications and ad hoc assumptions. First, we 
limit ourselves to the consideration of only one fiber direction and second, we 
consider a one- dimensional motion in the bulk of an infinite body. Here the 
motion is linearly polarized in a direction normal to the plane containing the 
direction of propagation and the direction of the fibers. We acknowledge that 
more complex anisotropies, geometries, and couplings arise in biomechanical 
applications. For instance, the mechanics of the aorta involves two families of 
parallel fibers, triaxial motions, and blood flow / arterial wall coupling. How- 
ever we argue that some major characteristics of biological soft tissues are 
encompassed in the choices of transverse isotropy, of infinite extent, and of a 
motion governed by an ordinary differential equation. Indeed the anisotropy 
due to the presence of one family of parallel fibers complicates the governing 
equations to an extent which is only marginally less than that due to the 
presence of two families of parallel fibers. Also, soft biological tissues are 
nearly incompressible and a (compressive) longitudinal wave is difficult to 
observe; it thus make sense to focus on transverse shear motions, which are 
useful in imaging technologies. Our third assumption is that the elastic strain 
energy is the sum of an isotropic part and an anisotropic part (called a rein- 
forcing model), in order to model an isotropic base material augmented by a 
uniaxial reinforcement in the fiber direction. Albeit strong, this constitutive 
assumption is now common and used by many authors (e.g. Triantafyllidis 
and Abeyaratne 1983, Qiu and Pence 1997, Merodio and Ogden 2002). Fi- 
nally, we assume that the solid is viscoelastic and here we assume not only 
Newtonian viscosity (proportional to the stretching tensor) but also fiber- 
oriented (anisotropic) viscosity. That latter assumption is strong but can be 
removed from our calculations by taking a constant to be zero. We believe 
that it might be useful in modeling the well-documented physiological effect 
of stretching training in sport medicine, which is that it affects the viscosity 
of tendon structures but not their elasticity (Taylor et al., 1990; Kubo et al., 
2002). 

We divided the article into the following sections. Section 2 presents the 
constitutive model and the derivation of the equation governing the rectilin- 
ear shear motion. As expected, this equation is nonlinear in the shear strain: 
it is a second-order partial differential equation, with cubic nonlinearity. To 



2 



initiate its resolution, we first look at the quasistatic experiment of recovery 
in Section 3. Then we have a first-order ordinary differential equation, and 
we find that it can lead to unusual behaviors when certain conditions (strong 
anisotropy, large angle between the shearing direction and the fibers) are 
met. The same is true of the case of creep, treated in Section 4. Basically, 
it turns out that the nonlinearity introduces ranges of material parameters 
and angles for which an expected behavior - say strain growth in creep - 
can be turned on its head - and lead to strain decay with time in creep, 
say. In the course of the investigation we develop synthetic tools of analy- 
sis which highlight the boundaries of these ranges. They also guide us for 
the resolution of the full dynamical equation of motion, which we tackle in 
Section 5 for traveling wave solutions. Again the solution may behave in an 
unexpected way, provided the anisotropy is strong enough and the fibers are 
in compression. Finally, Section 6 recaps the results and puts them into a 
wider context. 

2 Basic equations 

2.1 The viscoelastic anisotropic model 

We describe the motion of a body by a relation x = x(X, t), where x denotes 
the current coordinates of a point occupied by the particle of coordinates X 
in the reference configuration at the time t. 

We introduce F = <9x/<9X, the deformation gradient, and C = F T F, 
the right Cauchy-Green strain tensor. We focus on incompressible materials 
for which all admissible deformations must be isochoric, or equivalently, for 
which the relation det F = 1 must hold at all times. 

The body is reinforced with one family of parallel fibers. Our first as- 
sumption is that the unit vector a , giving the fiber direction in the reference 
configuration, is independent of X. The stretch along the fiber direction is 
Vao-Cao = V a ' a where a = Fao. 

We may now introduce the elastic part of our constitutive model. We 
consider the so-called standard reinforcing model, which is a quite simple 
generalization to anisotropy of the neo-Hookean model (Triantafyllidis and 
Abeyaratne, 1983; Qiu and Pence, 1997). For the standard reinforcing model, 
the strain-energy density is given by 



Here /i > is the infinitesimal shear modulus of the isotropic neo-Hookean 
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matrix, 7 > is the elastic anisotropy parameter, and the invariant J 4 
measures the squared stretch in the fiber direction. Mechanical tests show 
that the neo-Hookean strain energy function //(ii — 3)/2 fits uni-axial data 
rather well for arteries (Gundiah et al., 2006), while the anisotropic term 
70(^4 — l) 2 is adequate to describe a reinforced material which penalizes 
deformation in the fiber direction (Merodio and Ogden, 2003). 

The spatial velocity gradient L(X, t) associated with a motion is defined 
as L = grad v, where v = dx/dt is the velocity, and the stretching tensor 
D is defined as D = |(L + L T )- For incompressible materials, tr D = 
at all times. Newtonian viscous fluids possess a constitutive term in the 
form 2z/D, where v is a constant. For our special solid, we modulate the 
Newtonian viscosity with an anisotropic term, by replacing v with v[l + 
7i(/ 4 — 1)], where 71 > is the viscous anisotropy parameter. We show 
in the course of the paper that this simple choice of anisotropic viscosity 
captures the essential characteristics of attenuation in soft biological fibrous 
tissues. According to Baldwin et al. (2006), ultrasonic measurements of 
freshly excised myocardium show that "the attenuation coefficient was found 
to increase as a function of frequency in an approximately linear manner 
and to increase monotonically as a function of angle of insonification from a 
minimum perpendicular to a maximum parallel relative to the direction of 
the myofibers." 

We are now ready to give the complete Cauchy stress tensor of our vis- 
coelastic, transversally isotropic material as 

T = —pi + /i[B + 7o(/ 4 - l)a <g> a] + 2z/[l + 7l (/ 4 - 1)]D, (2.2) 

where the p is the yet indeterminate Lagrange multiplier introduced by the 
incompressibility constraint, and B = FF T is the left Cauchy-Green tensor. 

2.2 Shear motion 

We take a fixed, orthonormal triad of vectors (i, j, k), and call X, Y, Z the 
reference coordinates; hence X = Xi + Y j + Z\t. The triad is such that the 
unit vector in the fiber direction lies in the XY plane; hence, 

a = cos6>i + sin#j, (2.3) 

(say) where 6 G [0,7r] is the angle between the X-axis and the fibers. 
We then consider the rectilinear shearing motion, 

x = X + u(Y,t), y = Y, z = Z, (2.4) 
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where the anti-plane displacement u is real and finite. Then the components 
of the gradient of deformation F and of its inverse are given by 
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(2.5) 



where U = du/dY is the amount of shear. The left and right Cauchy-Green 
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(2.6) 



respectively, from which the expressions of the invariants I\ and J4 follow, 

I 1= 3 + U 2 , h = l + Usm26 + U 2 sm 2 6. (2.7) 

Figure shows the variations of I4 with 9 for several values of U between 
and 1. When I4 > 1 the fibers are in extension, and when J4 < 1 they are 
in compression; the figure shows that this latter behavior occurs in a smaller 
and smaller angular range, but is more and more pronounced, as the amount 
of shear is increased. Conversely, Figure [lb shows the variations of I 4 with U 
for several values of 9; when < 9 < tt/2, the fibers are always in extension 
and when 7r — tan -1 (2) = 2.034 < 9 < 71, they are always in compression 
for < U < 1. We refer to the paper by Qiu and Pence (1997) for similar 
figures and closely related discussions. 

In the deformed configuration, we find that a = (cos 9 + U sin 9)i + sin 9]. 
The remaining tensors required to compute the Cauchy stress tensor (12. 2p 
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(2.9) 



so that the non-zero components of T are T 33 = — p + /j, and 

Tn = ~P + //(l + U 2 ) + fxj {I 4 - l)(cos9 + Usm9) 2 , 
T22 = ~P + V + ^lo{h ~ 1) sin 2 9, 

Ti 2 = fxU + vj (h - l)(cos0 + Usm9) sin# + v[l + 71 (J 4 - l))U t . (2.10) 
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Figure 1: Variations of the squared stretch in the fiber direction: (a) with 
the angle and (b) with the shear. When J4 > 1, the fibers are in extension; 
when J4 < 1, they are in compression. 



Now the equations of motion div T = px tt reduce to the two scalar equa- 
tions: — p x + Ti2, y = putt and — p y + T 2 2, y = pu u - Differentiating the former 
with respect to y and the latter with respect to x, and eliminating p xy , we 
arrive at a single governing equation for the rectilinear shear motion, 

p U tt = nU yy + p-y sin 2 9 [U (2 cos 9 + U sin 9) (cos 9 + U sin 0)] 

+ vU tyy + i/7x sin0 [C/£/ t (2cos0 + t7sin0)] w . (2.11) 

Using the scalings t = ptjv and y = y/L (where L is a characteristic length 
to be specified later on a case- by-case basis), we write this equation in di- 
mensionless form as 

eU R = Uy y + 70 sin 2 9 [U{2 cos9 + U sin 9) (cos 9 + U sin 9)] 

+ U m + 71 sin 9 [UU t (2 cos 9 + U sin 9)] ~~ , (2. 12) 

where e = ppL 2 /v 2 . This is the main equation of our study. For convenience 
we drop the tildes in the remainder of the paper. We also introduce the 
following functions, 

f(j ,U,9) = l + 7 O sin 2 0(2cos0 + [/sin0)(cos0 + [/sin0), 
^( 7 i,C/,0) = l + 7 1 [/sin0(2cos0 + f/sin0), (2.13) 

so that (12.121) is now 

eU tt = [Uf(K, U, 9) + Utgbu U, 9)} vv . (2.14) 
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3 Nonlinear anisotropic recovery 



Our first investigation is placed in the quasistatic approximation, where we 
study the influence of elastic anisotropy and viscous anisotropy on the classic 
experiment of viscous recovery. We imagine that the material is sheared and 
that at t — the shear stress is removed: T 12 (0) = 0. Here the characteristic 
length L is the displacement at t — from which the material will relax to 
the unstressed state. 

In the quasistatic case, we neglect the inertia term of (I2.12j) and may 
thus integrate it twice to give the following first-order ordinary differential 
equation, 



Here we take the constants of integration to be zero, according to the context 
of recovery, as explained above. We then solve the equation as 



where the constant is computed so that U(0) = 1. 

When 9 = 0, the fibers are not active with respect to the deformation and 
we recover the classical result of isotropic viscoelastic recovery: U(t) = e~*. 

When 9 = tt/2, the anisotropic effects are at their strongest. In that case 
the integral above has a compact expression and we find: 



We now take 71 = (no anisotropic viscosity) and 70 = 1, 5, 100 (recall 
that the fibers are inextensible in the limit 70 — > 00). Figure [2|i shows that 
as the anisotropic effect becomes more pronounced, the recovery is quicker; 
in other words, the influence of elasticity becomes stronger as 70 increases. 
Then we fix 70 at 5 for instance, and look at the role played by the anisotropic 
viscosity, by taking in turn 71/70 = 0.5, 1.5, 2.5. We find in Figure [2b that, 
as expected, the viscous recovery is slower as 71 increases. 

When 9 7^ 0, 9 7^ 7r/2, other behaviors arise, which call for a detailed 
analysis. In particular, the exponential, or near-exponential, decay toward 
zero as t — > 00 is not necessarily ensured, especially when the anisotropic 
effects are strong and the fibers are oriented at a large angle 9 > tt/2. Clearly, 
Ut = when / = 0, according to (13.21) . Also, U t < when / and g are of 
the same sign and U t > when / and g are of opposite signs. These two 
functions are quadratic in U. If they have no real roots in U, then they are 



Uf( lo ,U,9) + U t g( ll ,U,9) = 0. 



(3.1) 




(3.2) 




(3.3) 
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Figure 2: Time recovery function when 9 = 7r/2: (a) 71 = and 70 = 1, 5, 
100; (b) 70 = 5 and 71 = 0.5, 1.5, 2.5. The recovery function for an isotropic 
solid is also plotted (dotted curve). 



both of the positive sign, and U t < (this is clearly the case in the region 
< 9 < 7r/2.) If they have real roots, then they may change sign and U 
might be an increasing function of t. This happens for / and for g when 
7r/2 < 9 < 7r and 

4 1 . . 

7o > — 2a . 2a , 7i > — 2A' 3 " 4 
cos 2 9 sin 9 cos 2 9 

respectively. In Figure the region C corresponds to the first inequality, 
where the delimiting curve has a vertical asymptote at 9 = 7r/2, a vertical 
asymptote at 9 = it, and a minimum at 9 = 7r/4, 70 = 16; we recall that 
Qiu and Pence (1997) showed that when 70 > 16, "simple shear at certain 
fiber orientations involves negative shear stress in the shearing direction for 
certain positive shears." The regions B and C correspond to the second 
inequality, where the delimiting curve has a vertical asymptote at 9 = tt/2 
and an horizontal asymptote at 71 = 1. 



3.1 Weak anisotropy 

First, we take both 70 and 71 in region A. This is the simplest case because / 
and g are then both positive and so U t is always negative (damped recovery). 
We took several representative examples in this region (say 9 = ir/4, 70 = 20, 
71 = 1) and checked, through integration and implicit plotting, that the 
graphs are indeed of the same nature as those in Figure EJ 
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Figure 3: Recovery: regions where the sign of U t may change. 
3.2 Strong elastic anisotropy 

Second, we take 71 in region A, by fixing it at 71 = 1, say. In that region, 
g > always and thus the sign of U t is the opposite of the sign of /. Then 
we take 70 = 20, which is above the minimum of region C. In Figure HI we 
plotted the locus for the values of U as functions of 9 such that /(20, U, 9) = 0. 
Outside the resulting oval shape, / > 0, and inside, / < 0. We also plotted 
the line U — 1, which intersects the oval at 6> min = 2.136 and # max = 2.221. 
Recall that U(0) = 1. 

When 9 > S max , U(t) starts at 1, and decreases because / > so that 
Ut < 0; as U decreases toward 0, U t tends to zero according to (]3.2jl i . but 
takes an infinite time to do so, according to (I3.2IW : hence U = is a horizontal 
asymptote and the recovery is "classical" , see plot (i) in Figure HI traced at 
6 = 2.4 (notice however that the recovery is not exponential because the 
second derivative of U clearly changes sign as t increases, in contrast with 
e - *, traced in dotted line.) 

When ^min < 9 < max , U(t) starts at 1 and then grows until it hits the 
upper side of the oval, taking an infinite time to do so; then this upper bound 
gives a horizontal asymptotic value, above the initial value, see plot (ii) in 
Figure IH traced at 9 = 2.2. 



9 



Figure 4: Types of time recovery functions for 7 = 20, 7i = 1 (strong elastic 
anisotropy). The amount of shear U starts at 1 for t = 0. Outside the oval 
shape, Ut < and U decreases as in (i) and (iv): decay toward zero; and 
in (in): decay toward a value > 0. Inside the oval shape, Ut > and U 
increases as in (ii): growth toward a value > 1. The recovery function e~* 
for an isotropic solid is also plotted (dotted curves). 

When 9 m < 9 < 6 min , where 9 m = 2.124 is the angle at which the oval 
plot has a vertical tangent, U(t) starts at 1 and then decreases until it hits 
the upper side of the oval, below 1 but above 0; then this lower bound gives 
a horizontal asymptotic value, above zero, see plot (iii) in Figure HI traced at 
9 = 2.125. 

Finally, when 9 < 9 m , U(t) starts at 1 and then decreases until zero; then 
this lower bound gives zero as a horizontal asymptotic value, see plot (iv) in 
Figure HI traced at 9 = 2.05. Notice that the second derivative changes signs 
three times as t increases. 
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3.3 Strong viscous anisotropy 

Third, we take 70 outside the C region, by fixing it at 70 = 1, say. In that 
region, / > always and thus the sign of U t is the opposite of the sign of 
g. Then we allow 71 to be in region B, and thus allow g (and Ut) to change 
sign with increasing 8, by taking 71 = 3.0 say. In Figure El we plotted the 
values of U as functions of 8 such that g(3,U,8) = and obtained the thick- 
line shape. Outside the shape, g > 0, and inside, g < 0. We also plotted 
the horizontal line {7=1, which intersects the shape at 8 m i n = 2.356 and 
#max = 2.820, and the vertical line 8 = 9 m = 2.186, which is tangent to the 
shape. 

Now when 8 < 8 m or 8 > # max , U(t) starts at 1 and decreases until zero; 
as U — > 0, the denominator in the integral tends to zero, indicating that it 
takes an infinite time to do so; hence, zero is a horizontal asymptote in these 
cases. To draw Figure E^i), we took 8 = 3.0 and for Figure EJ^iv) , we took 
8 = 2.1; both graphs show a somewhat classical decay with time. 

However, when 8 min < 8 < max , U(t) starts at 1 and then grows because 
U t > inside the thick line shape. Eventually U hits the upper face of the 
shape, where g = 0; then by ( 13. ip . either Uf = 0, or U t — > 00. Clearly, 
the first possibility is excluded because {7^0 when it is larger than 1, and 
/ 7^ when 70 is outside the C region. It follows that U grows and hits the 
upper face of the shape with a vertical asymptote after a finite time (and 
then stops because it cannot increase further since U t < outside the shape, 
it cannot remain constant since Ut 7^ on the shape, and it cannot decrease 
since U t > inside the shape.) Figure EJ^ii) shows such behavior for U(t), 
traced at 8 = 2.4. 

Finally, when 9 m < 8 < 8 min , U(t) decays from 1 until it hits the shape 
from above after a finite time; see Figure E^iii), traced at 8 = 2.2. Notice 
how quickly the final value is reached, compared to the isotropic exponential 
recovery. 

3.4 Strong elastic and viscous anisotropics 

In the case where both 70 and 71 are in the region C, any combination 
and overlaps of the thick curves presented in Figures H] and [5] may arise. The 
tools presented in the two previous subsections are easily transposed to those 
possibilities. A special situation arises when the locus of / = intersects the 
locus of g = 0; then, the numerator and the denominator in (13.21) may have 
a common factor so that the integrand simplifies and a regular behavior may 
appear. This situation is however too special to warrant further investigation, 
and we do not pursue this line of enquiry. 
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Figure 5: Types of time recovery functions for 70 = 1, 71 = 2 (strong viscous 
anisotropy). The amount of shear U starts at 1 for t = 0. Outside the thick- 
line shape, U t < and U decreases toward zero as in (i) and (iv). Inside 
the thick-line shape, U t > and U increases rapidly as in (ii), until it ceases 
to exist. There is also an angular region 9 m < 6 < # min where U decreases 
rapidly until it ceases to exist, see (Hi). The recovery function e - ' for an 
isotropic solid is also plotted (dotted curves). 
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4 Nonlinear anisotropic creep 



Our second investigation is again placed in the quasistatic approximation, 
where we now study the influence of elastic anisotropy and viscous anisotropy 
on the classic experiment of viscous creep. As the resulting analysis is similar 
to that conducted for recovery, we simply outline the main results. 

We imagine that the material is sheared and that the shear stress is 
maintained: Ti2(oo) 7^ 0. Here the characteristic length L is an asymptotic 
value of the displacement. We neglect the inertial term of ( 12.12ft and integrate 
it twice to give the following ordinary differential equation, 

Uf( l0 , U, 9) + U t g( lu U, 9) = const, (4.1) 

where we took the constant of the first integration to be zero, and the con- 
stant of the second integration to correspond to the applied (constant) shear 
stress, as is usual in the creep problem. More specifically, this constant is 
taken so that U(oo) = 1, and so is equal to /(70, 1,0); it follows that the 
equation above can be written as 

h( l0 , U, 9){U - 1) + 2(71, U, 6)U t = 0, (4.2) 

where h is defined by 

% , U, 9) = [l7/( 7 o, U, 9) - /( 7o , 1, 9)}/{U - 1) 

= 1 + 7 O sin 2 0[l + cos 2 9 + (U + 1) sin0(C/sin0 + 3cos0)]. (4.3) 
We then solve the equation as 



I 



( c/-i)M7».c,«) dC/ = - t + const " (44) 

where the constant is computed so that U(0) = 0. Hence the equations 
governing creep are almost identical to those governing recovery, with the 
difference that / is now replaced by h. 

Here we are mostly concerned with the question of how, if at all, a state of 
shear can be reached such that once removed, the unusual recovery behaviors 
of the previous section emerge. Thus we concentrate on strong anisotropic 
effects, with emphasis on strong elastic anisotropy (where the new function h 
is involved). We traced the regions where g and h, and thus U t , may change 
signs and found that the resulting graph is similar to that of Figure [3j with 
the main difference that the minimum of region C is now located at 9 = 37r/4 
and 70 = 4. Thus unusual behavior in creep may occur at much lower levels 
of elastic anisotropy than in recovery (where the minimum is at 70 = 16). 
We recall that Qiu and Pence (1997) show that when 70 > 4, "simple shear 
at certain fiber orientations involves a nonmonotonic relation between the 
shear stress in the shearing direction and the amount of shear." 
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4.1 Strong elastic anisotropy 

We begin with the case where h plays a major role, that is when 70 is greater 
than 4. For the purpose of direct comparison with the recovery problem, we 
take 70 = 20 and 71 = 1, as in Section 3.2. Figure [6] displays the curve where 
h(20, U, 9) = 0. Outside the thick-line curve, U t > 0, and inside, U t < 0. The 
curve intersects the line U = twice, at 9 min = 2.136 and at 6* max = 2.221. 
These are the values at which / = intersects U = 1 in Section 3.2 (see 
thin-line shape), because by (14. 3ft . h(20, 0, 9 min ) = /(20, 1, 9 min ) = 0, and 
similarly h(20, 0, 9 max ) = /(20, l,0 max ) = 0. We also display the vertical 
lines 9 = 9m = 2.664, where h = intersects U — 1, and 9 = 9 m = 2.042, 
where /i = has a vertical tangent. Recall that for creep, U(0) = 0. 

When 9 > 9m, U(t) starts at and grows toward 1; then U t tends to 
zero according to ( 14.21) . but takes an infinite time to do so; hence U = 1 is 
a horizontal asymptote and the creep is "classical" , see plot (i) in Figure El 
traced at 9 = 2.7 (the exponential creep function of isotropic visco-elasticity 
(1 — e _i ) is shown in dotted line.) 

When # max < 9 < 9 M or when 9 m < 9 < 9 min , U(t) starts at and then 
grows until it hits the oval shape, taking an infinite time to do so; then this 
upper bound gives a horizontal asymptotic value, below 1, see plot (ii) in 
Figure El traced at 9 = 2.5, and plot (iv), traced at 9 = 2.05. 

When 9 min < 9 < 9 max , U(t) starts at inside the oval shape and thus 
it decreases until it hits the lower side of the shape, taking an infinite time 
to do so; then this lower bound gives a horizontal asymptotic value, below 0, 
see plot (iii) in Figure El traced at 9 = 2.17. 

Finally, when 9 < 9 m , U(t) can again grow toward 1, see plot (v) in Figure 
El traced at 9 = 2.0. Notice however that the concavity of the curve changes 
as t increases. 

4.2 Strong viscous anisotropy 

Here we remark that the function governing the strength of the viscous 
anisotropy, namely g, is the same for creep as it is for recovery. Thus, the 
region where U t might change sign because of strong viscous anistropy is 
the region B of Figure [2J Also, the locus of points where g = is typically 
displayed by the thick-line shape of Figure [5] and because g(ji, 0, 9) = 1 > 
always, this curve never crosses the abscissa U = 0. It follows that there is 
only one situation where viscous anisotropy leads to anomalous creep, when 
#mm < 9 < 9 max ; then, U(t) starts at zero, and grows toward the thick-line 
shape, which it reaches after a finite time with a vertical asymptote. 
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Figure 6: Types of time creep functions for 70 = 20, 71 = 1 (strong elastic 
anisotropy). The amount of shear U starts at for t = 0. Outside the thick- 
line shape, U t > and U increases as in (i) and (v): growth toward 1, and as 
in (ii) and (iv): growth toward a value below 1. Inside the thick-line shape, 
U decreases as in (iii): decay toward a negative value. The creep function 
1 — e _i for an isotropic solid is also plotted (dotted curves). 
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4.3 Pre-stretch and nonlinear anisotropic creep 



Here we show how anomalous creep can be avoided (amplified) by stretch- 
ing (compressing) the solid prior to the shear. Hence, instead of ( 12.41) . we 
consider the motion 



x 



A~2X + Xu(Y, t), y = XY, 



A"27 



(4.5) 



The following decomposition of the associated deformation gradient shows 
that the solid is stretched by a ratio A in the Y direction, 
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(4.6) 



(Note that F 2 Fi ^ FiF 2 .) The kinematic quantities of Section 2.2 are 
modified accordingly. In particular, 



A" 1 cos z 6 + X z sin^ 9 + U\2 sin 26 + U' Z X Z sin 



2\2 



(4.7) 



The end result is that the differential equation governing creep is changed 
from (T£2D to 

h\ l0 , U, 9)(U - 1) + 2 A ( 7 i, U, 9)U t = 0, (4.8) 
where h x and g x are defined by 

h x ( l0 ,U,6) = A 2 {1 +7 O sin 2 ^[2A 2 sin 2 + 3A" 1 cos 2 0- 1 

+ (U + l)sin^([/sin^ + 3cos^)]}, 

g x ( lu U,6) = l+ ll {\~ 1 cos 2 6 + \ 2 sm 2 6-l + U\hm28 + U 2 \ 2 sm 2 6). 

(4.9) 

Figure [7] shows the loci of h x = in the case of a strong elastic anisotropy 
(70 = 30, 71 = 1), for several values of A. The figure clearly shows that the 
pre-stretch A can be used to control the shape of these curves: if the solid 
is put in compression first, and sheared for creep next, then the region of 
potential anomalous creep is increased; if it is put under tension, then the 
area of the region rapidly decreases and eventually dissapears altogether. 



5 Nonlinear traveling waves 

So far we have looked at how the presence of elastic and viscous fibers affects 
some quasi-static processes. Typically, creep and recovery connect one state 
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Figure 7: Effect of pre-stretch on time creep functions for 70 = 30, 71 = 1 
(strong elastic anisotropy). The amount of shear U starts at for t = 0. 
Inside the thick-line shapes, U decreases; this situation may arise when the 
solid is not pre-stretched (A = 1) or when it is compressed (A = 0.8). Outside 
the thick-line shapes, Ut > and U increases; this situation may arise when 
the solid is in extension (A = 2, A = 2.7). 
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of constant shear (initial) to another (final). Now we examine another class 
of solutions connecting two constant states of shear, this time dynamically, 
by looking for traveling wave (kink) solutions. 

The mathematical theory of one-dimensional transverse traveling waves 
in isotropic viscoelastic materials with a Kelvin- Voigt type of constitutive 
equation is well grounded, see for example Nishihara (1995) for a clear and 
complete mathematical approach, or Jordan and Puri (2005) for a specific 
and explicit example. A traveling wave is a solution to the equations of 
motion in the form 

U(Y,t) = U(£), £ = Y-ct, (5.1) 
where c is the constant speed; also, U is such that 

lim U(0 = U L , KmU(0=U R , (5.2) 

where Ul and Ur are distinct constants. In what follows, we focus on the case 
where Ul = 0, Ur = 1. This case is general up to a rigid translation. Here 
we take the displacement corresponding to Ur as the characteristic length L. 
Substituting (15.11) into (13. 2p . we obtain 

ec 2 U" = (Uf - dJ'gf , (5.3) 

and then by integration, 

cU'g = (f - ec 2 ) U + const. (5.4) 

By the requirement Ul = 0, the constant must be zero. By the requirement 
Ur = 1, we have 

/( 7o , 1,9) = ec 2 . (5.5) 

This equation prompts three remarks. 

First, we must ensure that /(70, 1, 9) > 0. Recall that according to (12.131) . 

/(7o,l,0) = 1 + 70 sin 2 9(2 cos 9 + sin 6) (cos 9 + sin 9), (5.6) 

and so, 

5/(70, 1, 9) Id9 = 70 sin 9(4 cos 3 9 + 9 sin 9 cos 2 9-3 sin 3 9). (5.7) 

In Figure [HK we plotted the variations of [/(70, 1,9) — l]/7o with 9, as well 
as those of its derivative with respect to 9 (scaled to 1/8). Clearly, the 
function ( 15.61) . viewed as a function of 9, has an absolute minimum and an 
absolute maximum. The minimum is at 9, say, such that tan# is that root 
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Figure 8: (a) Variations with 9 of [/(70, 1,0) — l]/jo and of its derivative, 
showing an absolute minimum at — 2.1777. (b) Variations of —3/ tan# — 1 
with 9, crossing the abscissa line at #0 = 1.8926. 

of the cubic 4 + 9x — 3x 3 = corresponding to n/2 < 9 < n; numerically, 
9 = 2.1777. Then, solving /( 7o , 1, 9) = for 7o , we find that /(70, 1,9) > 
when < 70 < 70 = 18.490; and that when 7 > 70, there appears a range 
for 9 where /(7o,l,#) > is not insured. Placing ourselves outside that 
possibility, we deduce from (15. 5ft that for a given 70 and a given 9, the wave 

travels with speed 

c = ±V/(7o,l,0)/e. (5-8) 

This is of course expressed in the dimensionless variables of length/L and 
timex fi/u. Turning back, if required, to physical variables, we would find 
that the wave travels with the dimensional speed \J /i/(7o, 1, 0)1 v. 

The second remark is that according to ( 15. 5ft and ( 15.61) . the wave (when 
it exists) travels with maximum speed at the angle 9, say, such that tan# is 
that root of the cubic 4 + 9x — 3x 3 = corresponding to < 9 < tc/2; numer- 
ically, 9 = 1.0910. Hence the directions of extremal speeds of propagation 
are always the same, whatever the values of the constitutive parameters /i, 
70, and 71 are. This observation indicates the way for an acoustic determi- 
nation of the fiber orientation: if an experimental measurement of the shear 
wave speed can be made in every direction of a fiber-reinforced viscoelastic 
nonlinear material, then the fibers are at an angle 9 from the direction of the 
slowest wave and at an angle 9 from the direction of the fastest wave. We re- 
call that for waves in an isotropic deformed neo-Hookean material, Ericksen 
(1953) found that the fastest waves propagate along the direction of greatest 
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initial stretch. 

The third remark is that when (15.51) holds, then 

/(70, U, 9) - ec 2 = 7 U(U - 1) sin 3 9 [U sin 9 + 3 cos 9 + sin 9} . (5.9) 

Then the separation of variables, followed by integration of the first order 
differential equation (15.41) . leads to 

f : 9(lfi,U,0) = to cong ( 1Q) 

J U(U -1) sin 3 9[Usin9 + 3 cos9 + sin9] c v ' 

where the constant of integration is arbitrary; without loss of generality, we 
take it to be such that U(0) = 1/2. 

Clearly, critical issues arise when either the numerator or the denominator 
change signs (because then U' changes sign and it might not be possible 
to find a solution satisfying the requirements ( 15. 2ft ). We may take care of 
the numerator's sign by considering elastic anisotropy only (70 7^ 0) and 
discarding viscous anisotropy (71 = 0); then g — 1. For the denominator 
however, we note that U sin 9 + 3 cos 9 + sin 9 can change sign for certain 
ranges of U and 9. Figure [Hb shows the curve U = —3/ tan 6* — 1; on its left 
side, the denominator is positive, on its right side, it is negative. Accordingly, 
the wave connects to 1, see in Figure[9|i, or is unable to do so, see Figure[9b. 
In that latter case, the wave front grows toward an asymptotic value which 
is less than 1; a second solution exist (dotted curve) with 1 as an asymptotic 
value, but in the £ — > —00 direction. 

As a final remark, we note that when 71 is large enough to allow for the 
possibility that g = (strong viscous anisotropy), then a "singular barrier" 
arises, see Pettet et al. (2000). 

6 Discussion 

In the course of this investigation on nonlinear anisotropic creep, recovery, 
and waves for fiber-reinforced nonlinear elastic materials, we unearthed some 
complex mechanical responses. For some range of the constitutive parameters 
and for some angle ranges of the fiber arrangement, we saw that unusual and 
possibly aberrant behaviors can emerge. 

From a mathematical point of view, we gave a detailed explanation of 
the reasons for these behaviors, by linking them to the singularities of the 
determining equations for the amount of shear. 

From the mechanical point of view, we pointed out that non-standard 
behaviors always occur when the angle between the fiber family and the 
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Figure 9: Traveling wave solution for anisotropic elasticity (71 = 0); (a) at 
9 = 1.8, (b) at 9 = 2.1. 

direction of shear is such that the fibers are compressed, see Figure [H It 
has been widely demonstrated that several types of instabilities may develop 
in the case of fiber contraction, see the detailed studies by Triantafyllidis 
and Abeyaratne (1983), Qiu and Pence (1997), Merodio and Ogden (2002, 
2003, 2005ab), or Fu and Freidin (2004). For example, Merodio et al. (2007) 
recently investigated a non-homogeneous rectilinear shear static deformation 
for the standard reinforcing model (12.1 1) and found non-regular solutions (that 
is, deformations characterized by a discontinuous amount of shear) in fiber- 
contracted materials. 

From a numerical point of view, we recall a simple model, together with 
a simple class of solutions, allows a step-by-step control of the simulations. 
It would indeed be hard to detect non-standard behaviors by relying solely 
on a complex numerical finite element method, and omitting to conduct a 
simple analytical methodology such as the one presented in this paper. For 
example, Holzapfel and Gasser (2001) present a detailed computational study 
of some viscoelastic fiber-reinforced nonlinear materials, but use values for 
the material parameters and for the angles which place their simulations 
outside the problematic ranges. Other studies are placed in the framework 
of linear models (even for polymeric materials, see Liu et al., 2006), which 
fail to capture non-standard behaviors. 

From an experimental point of view, our results suggest some simple yet 
revealing protocols. In particular, it would be most valuable to investigate the 
existence and the persistence of asymptotic residual shear strains, sustained 
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after the shear stress is removed, at levels not only below the value at initial 
time but also above (as in Section 3). So far we have only identified reports 
of experimental results concerned with elastomeric materials reinforced with 
inextensible fibers (and therefore with a ratio between the shear modulus of 
the bulk matrix and that of the fibers of several orders of magnitude), or 
concerned with moderate angles between the direction of shear and the fiber 
direction. 

From a biomechanical point of view, the results have meaningful impli- 
cations for biological soft tissues. First, the model captures adequately the 
elastic and the viscous anisotropies of biological materials (Baldwin et al. 
2006, Taylor et al. 1990). Second, although anomalous creep behaviors 
might preclude anomalous recovery behaviors, it is still useful to study the 
latter, because they might nonetheless arise in vivo following a stress-driven 
fiber orientation remodeling (Hariton et al., 2006). Third, the effect of the 
prestretch on non-standard behaviors is significant theoretically (Section 4.3) 
as well as practically (in vivo experiments show that large static prestretches 
of tendons reduces the risk of unexpected behaviors, see Kubo et al. (2002)). 
Finally, the results of the traveling wave study (Section 5) may eventually 
lead to an acoustic (elastographic) determination of the fiber angle in soft 
tissues, through an efficient, simple, and non-invasive investigation. 

Obviously, our results must be improved and several directions are pos- 
sible. Hence, two families of fibers have to be considered to give a better 
comparison with in vivo results for soft tissues. Also, the more realistic 
models of fiber reinforcements (such as the one proposed by Horgan and 
Saccomandi (2005) and by Gasser et al. (2006)) must be incorporated in the 
present study, to identify with a greater precision the range of parameters 
for which strange behaviors may occur. 
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